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Abstract 

Traditional numerical discretizations of conservative systems generically yield an artificial sec¬ 
ular drift of any nonlinear invariants. In this work we present an explicit nontraditional algorithm 
that exactly conserves these invariants. We illustrate the general method by applying it to the three- 
wave truncation of the Euler equations, the Lotka-Volterra predator-prey model, and the Kepler 
problem. This method is discussed in the context of symplectic (phase space conserving) integra¬ 
tion methods as well as nonsymplectic conservative methods. We comment on the application of our 
method to general conservative systems. 
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I. Introduction 


For many years now symplectic integrators have been the subject of much productive study. 
(See Channell and Scovel^^] for an overview; see also the recent book by Sanz-Serna and 
Calvot^l) There are variety of Hamiltonian systems for which symplectic methods have 
proven extremely useful, if not essential; but these methods do not constitute the last word on 
integration techniques. As Ge and Marsden show,[^] exact energy conservation is, in general, 
not possible with a symplectic method. Since the energy error is typically not secular but 
rather oscillatory, it is commonly believed that exact energy conservation is not as important 
a benefit as preserving the phase space structure. 

Concerning the numerical preservation of more general constants of motion, less is known. 
Based on the work of Cooper,W Sanz-Sernat^’®] has shown that a restricted class of quadratic 
invariants will be conserved by certain symplectic Runge-Kutta schemes. For the Runge- 
Kutta methods studied by Cooper,W conservation of quadratic invariants necessarily requires 
that the method be implicit. One method of ensuring the preservation of any constant of 
motion is to use the constant to reduce the number of equations that must be solved. If 
the constants are in involution, then an entire degree of freedom (one coordinate and one 
momenta) can be removed from the dynamics for each such constant. This is seldom practical 
since the relationship between the constants of motion and a given dynamical variable may 
well be noninvertible (see the discussion in Gear!®]). The net result is that the reduced 
equations tend to be more complicated than the original system (hence the “force” terms are 
more expensive to compute); thus, in a system with a large number of degrees of freedom, 
little advantage is gained. Furthermore, if the constants of motion are not in involution, the 
system obtained by eliminating these invariants will be noncanonical,[^’®] resulting in even 
greater complexity. 

It may be that the system of interest is most naturally described by variables that give 
rise to a noncanonical Hamiltonian structure. For noncanonical systems, Ge and Marsdent^l 
have provided a general construction for integrators that preserve both momentum maps and 
the structure of the Poisson manifold. Channell and Scovel^^l have shown how to implement 
these algorithms without the need of coordinatizing the configuration space group. This 
notwithstanding, they report that, with the exception of certain special (albeit important) 
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forms of the Hamiltonian, such methods tend to be computationally expensive. 

There is a further class of dynamical systems that is of interest, namely those systems 
that are not Hamiltonian (canonical or otherwise) but still possess constants of motion. 
Prime examples of such systems are transport equations, such as the Boltzmann equation. 
Further examples are truncations of the Fourier-transformed Euler fluid equations. The 
untruncated equations constitute an inflnite-dimensional Hamiltonian held theory, however, 
when the number of Fourier modes is reduced to a flnite set, the Hamiltonian structure is 
typically lost even though energy and enstrophy are still conserved. One might (justiflably) 
argue that such a truncation is not appropriate, but presently there is no practical alternative. 
(It is very much an open question as to the overall effect on the dynamics of such truncations.) 
Given that these systems are not Hamiltonian, symplectic methods are of no use, while the 
preservation of constants of motion is still of great interest. 

A variety of methods for enforcing conservation of general invariants has been proposed. 
Baylis and Isaacsont^^’^^l have proposed a two stage algorithm where the approximate solu¬ 
tion, obtained by standard methods in the first stage, is projected onto the constraint surface 
defined by the invariants in the second stage. Brasey and Hairert^^l have proposed a “half¬ 
explicit” method where the projection [via a Lagrange multiplier) and integration stages are 
merged together. LaBudde and Greenspant^^’^^l have developed an algorithm for central 
force problems that conserves both energy and angular momentum. Gear^®’^®] advocates an 
approach that amounts to an embedding of the original system into a higher dimensional 
space, yielding a set of differential-algebraic equations, the solution of which coincides with 
the solution of the original equations and preserves the invariants. 

Our purpose in this paper is to present another approach to the development of exactly 
conservative algorithms. We begin with a simple model problem, that is of interest in both 
fluid mechanics and plasma physics, which possess two quadratic invariants. We develop 
integrators for this system that are explicit and exactly conserve both invariants. We further 
illustrate our method by applying it to the Lotka-Volterra predator-prey model and to the 
Helper problem. 
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II. A Model Problem 


Our original interest in the issues of exact preservation of constants of motion arose in 
the stndy of two-dimensional inviscid flnid tnrbnlence. As an illnstration, consider the 
“three-wave” problem obtained by restricting the Fonrier-transformed Enler eqnations to 
three modest^®^^®]; 


d^l)^ 

dt 

d'lpp 

dt 

dt 


MriI^piI^q — Sp-{'ip) , 

(la) 

'^Q '^K — 5 

(16) 

'^P = ’ 

(Ic) 


where 'll) = ('0^,'0p,'0g), K, P, and Q are the magnitndes of the Fonrier wavennmbers of 
the three modes and the mode conpling coefficients Mp and Mq satisfy 


Mp Mq — 0 , 


( 2 ) 


and 

K^Mj^PP^MpPQ^Mq = 0. (3) 

This system possess two invariants: the total energy 

^ ~ 2 + '*/’p + (4) 

and the total enstrophy 

z = l{K^<Pl + P^<Pl + Q^i>l)- ( 5 ) 

The constancy of these qnantities follows directly from properties of Sf,: 

= (6«) 

k 

= (6^) 

k 

where k ranges over the set {K,P,Q}. (These equations are isomorphic to Euler’s equations 
for the rigid body, in which case the second invariant is the total angular momentum.) 
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When (1) is integrated numerically using standard methods neither E nor Z are exactly 
conserved. This behavior is made apparent by applying Euler’s method with a time step r: 

+ + ke{K,P,Q}. (7) 

The energy at the new time is 

E{tPT) = 

k 

k 

= E(t) + E^Y.Sl ( 8 ) 

k 

where we have used (6a) in the last step. Thus the total energy is always increasing. A 
similar calculation for the enstrophy gives 

Z(i + T) = Z(«) + iT=y]*^S^ (9) 

k 

which is likewise always increasing. For extremely long runs these results imply that a very 
small time step is required to keep the accumulated error down to a given level — clearly an 
undesirable situation. 

Many authors have noted that the lack of preservation of constants of motion potentially 
introduces significant unphysical effects and as such these errors are, in some sense, more 
important than those numerical errors that do not alter constants of motion. As de Frutos 
and Sanz-Sernat^^^ point out, one can think of the local error in a numerical integration as 
having two “components”: one which leads to unphysical changes in the constants of motion 
and another which does not. When these local errors accumulate over many time steps, the 
former component is significantly more harmful than the latter in that errors which lead to 
changes in the constants of motion affect the qualitative nature of the solution, whereas other 
errors only affect the quantitative results. (A similar observation regarding the accumulation 
of error in area-preserving maps has been made by Greene.In essence we are saying 
that nonconserving integrators have the potential to make “structural” errors in the solution 
— an observation which agrees well with one’s physical intuition. In the context of our 
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model problem the implication is clear: keeping the time step small enough to maintain a 
reasonable level of energy and enstrophy conservation means that we are likely to be using 
more computational resources to obtain a given accuracy in the solution than would otherwise 
be necessary with a conservative integrator. 

Although the three-wave problem is both integrable and Hamiltonian, our ultimate in¬ 
terest in this problem concerns the n-wave generalization of this system, which possesses 
both energy and enstrophy invariants but is not Hamiltonian; hence, we are lead to consider 
methods that do not rely on a particular geometrical structure. One might be tempted to 
enforce energy and enstrophy conservation by using these invariants to eliminate two modes 
from the dynamics. In this case the algebraic relations are simple enough to allow this, but 
there is a compelling physical argument against this approach. The modes are associated 
with different length scales; the choice of which modes to eliminate in favor of the invariants 
therefore has significant physical implications. Furthermore, the numerical error responsible 
for the lack of energy conservation could be imagined to contribute to a nonphysical energy 
and enstrophy transport. In any event, since each mode makes a positive contribution to 
both the energy and enstrophy, this scheme would eventually fail when the artificial growth 
in the invariants surpasses the initial contributions of the eliminated modes. Doubtless this 
scheme would exhibit nonphysical behavior long before this stage of failure was reached. 

III. Conservative Integrators for the Model 
Problem 

In light of the above discussion, an algorithm that exactly conserves energy and enstrophy is 
clearly desirable. As we have noted in Section I, a variety of implicit methods are known that 
preserve quadratic invariants. While implicit methods have noteworthy stability properties, 
they tend to be less computationally efficient than explicit methods since they typically require 
multiple evaluations of the “force” terms. Therefore, we turn our attention to the development 
of explicit conservative methods for our model problem. 

An elegant approach to this problem is found by borrowing from the ideas of backward 
error analysis.The essential idea is to construct a new system of equations that, under 
the conventional (nonconservative) integrator, yields a conservative numerical approximation 
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to the original equations, 
equations of the form 


To this end, consider the alternative problem described by three 


dt 


SM + ft- 


( 10 ) 


Our objective is to find an that guarantees exact energy and enstrophy conservation and 
that vanishes in the limit of small step size. The form of will depend on the integration 
algorithm. We begin by deriving for Euler’s method. We then construct a second-order 
predictor-corrector scheme. 


A. Euler’s Method 


As a “proof of principle” test we develop a conservative version of Euler’s method. While not 
particularly useful in practice, Euler’s method has the advantage that the algebra associated 
with constructing the conservative method is quite straightforward. 

Application of Euler’s method to the modified system yields 


+ = M^) + ^iSk + fk)- 


(11) 


The energy at the new time, 

Eit + r) = + 

k 

= E{t) + 2 [‘^^h'^k + ^^i^k + fk)^] ’ 

k 


will be conserved provided 


( 12 ) 


+ + = 0 - ( 13 ) 

k 

There is considerable freedom in satisfying (13). To guarantee that our approximate solution 
satisfies the original differential equation, it is necessary that vanish as r —0. That is, 
for small time steps, we must recover the original integration algorithm. One would prefer 
that ff, not introduce additional couplings into the differential equations. In light of these 
observations, let us try to satisfy (13) with the more restrictive condition that each term in 
the sum must independently vanish: 


24^fc + ^('5fc +4)^ = 0. 


(14) 
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There is an additional motivation for splitting (13) into three equations, namely that for 
satisfying (14), the enstrophy will also be conserved. This equation is easily solved, yielding 

^fk = + + + (15) 

where = a^(t,T) is so far an unknown sign. Evaluation of (15) at r = 0 implies 
that (T^(t,0) = sgn{'ip Upon substituting (15) into the Euler integrator, (11), we ob¬ 
tain the following time stepping rule: 

+ + (16) 

It is now clear that (T^(t,r) must in fact be the sign of -\- r). 

If '0^(t) ^ 0, then for sufficiently small r the sign can be expressed explicitly as = 
sgn{'tpf,{t)). In this limit, then vanishes, or equivalently, (16) reduces to Euler’s method: 

+ ^) = sgn(^fc(t ))+ 

-^k + ^^k- (17) 

In this case the new algorithm predicts values of 'ipf,{t + T) that are quite close to those given 
by Euler’s method — this is exactly what one would expect. The energy and enstrophy 
errors arising from (7) are the result of small (but nontrivial) errors in 'ipf^{t + T) that can be 
corrected by making a slight modification to the algorithm. 

However, if = 0, it is seen from (15) that has the nonzero limit —Sf, as r —?• 0. 
Consequently, (16) has a spurious fixed point at + t) = 0. Eortunately, a remedy for 

this problem has been developed. Equation (16) may be used up to and including the time 
step where 'ip{t -|- r) = 0. Likewise, a scheme that steps backwards in time from t 2 t may 
be used to integrate all the way back to t -|- r, where the two solutions must match. In this 
manner, we have developed an algorithm for the case where changes sign. 

Another potential problem with (16) is that the argument of the radical can become 
negative. In this case, we can rewrite the radical term as Xky where + 2 r 5^ 

is just the Euler approximation for a step size of 2r. The condition 'ip^Xk <0 tells us that 
Euler’s method is predicting a sign change of between t and t -|- 2t; hence, we are in 
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the vicinity of tjjj, = 0. To deal with this case, one adjusts the time step so that at the new 
time '0^ = 0 and then proceeds with the implicit algorithm described above. 

We give the name “Conservative Euler” (C-Euler) to this combined algorithm. In Eig- 
ure 1 we compare the numerical solutions of the three-wave problem obtained using the 
conventional Euler method with those obtained using C-Euler and with the exact solution. 
Eor these calculations K = Vs, P = S, Q = \/6, = 1, Mp = 1 and Mq = —2. The 

effect of the energy growth can be seen on the amplitudes computed by the Euler method. 

With the exception of the points where changes sign, C-Euler is an explicit algorithm. 
We will soon see that the gymnastics associated with the fixed point of the integrator just 
described are a consequence of the low order of the Euler method and that a fully explicit 
conservative integrator is possible. 



Figure 1: Solutions of the three-wave problem for the initial conditions = 
VP5, ipp = 0.0, and 'ipQ = vTTb computed using the conventional Euler 
(Euler), the conservative Euler (C-Euler) and the exact solution (Exact). 
A fixed time step of 0.02 was used. The unphysical energy growth in the 
conventional Euler algorithm leads to large errors in the amplitudes. 
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B. Predictor—Corrector 


In practice, one would prefer to use a scheme that both is of higher order than Euler’s method 
and has better stability properties. We now turn to a simple second-order predictor-corrector 
scheme, which we apply to our model problem (1): 


^k = ^k + ^Sk 

Mt + ^) = ^k + ^{Sk + Sk) , 


(18a) 

(186) 


where Sf, = Ski'ip)- As we will show, using a second-order method overcomes the fixed-point 
problem that we encountered with Euler’s method. 

The energy now evolves according to 


E{t + t) 


Y, + ri>k{Sk + Sk) 


E(t) + 

m + 







k 


7-2 / ^ \ 2 

+ {Sk + Sk) 


(19) 


where we have used the definition of ^pk and the properties of Sf, in the final step. A similar 
calculation gives 

Z{t + T) = Z(t) + ^'^e[s,-S,y. ( 20 ) 

k 

Again we see that the numerical method yields an ever increasing energy and enstrophy.^ 

To obtain a conservative version of this algorithm, we proceed as above by applying the 
predictor-corrector method to the modified equation of motion, (10), giving 


i^k = i^k + ^ i^k + fk) > 

+ '^) = + 2 + fk + ^k + h) ■ 


(21a) 

(216) 


^ One might be tempted to conclude that any conventional method will yield a positive-definite energy 
growth. While nonconservation is generic, the sign of the energy error is typically indefinite. For example, 
a second-order Runge-Kutta method gives oscillatory errors in energy and enstrophy, although on average 
both the energy and enstrophy grow. 
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As we commented above, the conservative algorithm makes only small corrections to the 
values of + This immediately brings to mind the underlying philosophy of the predic¬ 
tor-corrector algorithms; in fact, one might suspect that energy and enstrophy conservation 
can be achieved by modifying only the corrector part of the integrator. Since the predictor 
is merely an intermediate approximation, there is surely no need for it to be conservative. 
Thus we can replace (21) with the simpler prescription 


4 = 4 + (22a) 

+'^) = + 2 • (22^) 

As before, we determine g^, by demanding conservation of energy and enstrophy. The 
energy at t -|- r is given by 

E{t + T) = - + Ti>f, (^Sf, + + — (^Sf, + Sf, + 

= m + ^Y.[9k^k-rS,S, + ^- {S, + S, + g,f] , (23) 

k 

where the last step follows from the definition of the predictor and the properties of We 
see that energy will be conserved provided that 

E +1 (St + St + at)"] =0- (24) 

k 

Similarly, enstrophy will be conserved if 

-b - -b -b 5'^)^ =0. (25) 

k 

We can satisfy these conditions simultaneously if we can solve 

-b - (5'^ + 5'^ + 5'^)^ = 0 (26) 

for g^. Some straightforward algebra gives 
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where we choose = ±1 such that as r —> 0, vanishes. We consider the limit of small r 
in two cases. If '0^ is nonzero, then for small enough r, both and have the same sign 
and we can expand the radical to give 


T 


'9k — 


{^k + ^k) + ^k ^ + ^k) ] + 0{t^) ’ ( 28 ) 


leading us to choose = sgn(';/)^). Otherwise, if '0^ = 0, then pjf, = tS f^ and 5^ = <S'^ + 0(r), 
so that 


l9k = -^Sf^ + <^k\/^^ + 0{r^) 

=+ sgn(5J 5;- +0(r2). (29) 

In this case we take = sgn(S'^) = sgn('0^). In the previons case, we noted, for small r, 
that and have the same sign. Therefore, the choice = sgn('0^) will always provide 
the correct limiting behavior. 

Using the expression (29) for g^, in our modified predictor-corrector algorithm, (22), we 
obtain the following conservative integrator: 

^k = '>Pk^^Sk^ 

+ ^) = ^k \J'^k + ^{ASk + ^kSk) > 

where = sgn('0^). Unlike the C-Enler algorithm, this algorithm, which we call “conserva¬ 
tive predictor-corrector,” (C-PC), does not snffer from fixed points. It is still possible that 
the argnment of the radical can become negative; however, this merely indicates that the step 
size is too large. 

We now compare the numerical solutions of our model problem obtained with the con¬ 
ventional predictor-corrector method with those obtained from C-PC, (30). Our results are 
summarized in Figures 2-5. In Figure 2 we show computed with both methods as well 
as the exact solution. The errors in the two approximate solutions are displayed in Figure 3. 
The conservative predictor-corrector is seen to yield a more accurate solution. In Figure 4 
we plot AE = E{t) — F'(O) and AZ = Z{t) — Z{0) for both methods. 


(30a) 

(306) 


12 




Figure 2: Solutions of the three-wave problem for the initial conditions = vTTh, 
ipp = 0.0, and 'ipQ = vTTh computed using the predictor-corrector (PC), the 
conservative predictor-corrector (C-PC), and the exact solution (Exact). A fixed 
time step of 0.2 was used. 
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Figure 4: Change in energy and enstrophy for the conventional predictor-corrector 
method (PC) and the conservative predictor-corrector (C-PC). 


10 -3 


•S 10-4 
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10-3 
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Figure 5: Single-step error in mode K for the initial conditions = vTTS, 
ipp = 1.0, and 'ipQ = y/TPE for the conventional and conservative predictor-cor¬ 
rector methods. The results of fitting the error to a power law Ar" is shown, 
indicating that conservative algorithm is of second order, as expected. 
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In the limit of small step size, C-PC reduces to the conventional predictor-corrector. To 
illustrate this property we fit the error for a single step of mode K to a power law: 

Aij^ = AT^. (31) 

The results of this fit, shown for both methods in Figure 5, are consistent with our expectation 
that both the conventional and conservative predictor-corrector methods are second order. 

In constructing our conservative algorithms, we have essentially altered the manner in 
which truncation error enters the solution. Where this error has gone is an important question. 
It is unreasonable, of course, to expect that the truncation error has vanished. In fact, all 
of the truncation error is now lumped into the only place left — the phase of the solutions. 
Since our ultimate application is to fluid turbulence, the nature of this phase error could be 
of great importance. There are two general possibilities: either this error manifests itself as 
a global phase shift, with all three waves exhibiting the same phase error with respect to 
the exact solution, or each wave receives a different phase error, so that relative phase shifts 
begin to develop. Of the two possibilities, the first is of little consequence in a turbulence 
simulation, whereas the second could, arguably, be as bad (from a structural point of view) 
as the energy growth that we have sought to eliminate. 

Since our model problem is integrable, we can easily distinguish between these cases. In 
a plot where each of the axes is one of the dynamical variables, an integrable system yields a 
simple closed curve. For our case, such a plot is shown in Figure 6. The solid line is the orbit 
computed with the conservative integrator while the dots represent the solution obtained from 
the conventional predictor-corrector. Since the conservative solution yields a closed curve, 
we may conclude that the additional phase error introduced is global and thus the relative 
phases of the waves are not affected by our method. This supports the general observation 
made by de Frutos and Sanz-Serna regarding the nature of local truncation error in systems 
with invariants. 
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gure 6: Integration of the three-wave problem using a conve 
der predictor-corrector (solid line) and the conservative prec 
ots). The effect of the 4% energy gain by the conventional m 
sible. 


C. Generalizations 


The C-PC algorithm has two important straightforward generalizations: to n-waves and to 
complex The n-wave generalization is immediate — nowhere in our derivations of the 
conservative algorithms have we made use of the number of modes. Both C-Euler and C- 
PC can be applied to a system with an arbitrary number of modes, where the energy and 
enstrophy expressions are the appropriate generalizations of (4) and (5) respectively. 

The generalization to complex amplitudes proceeds as follows. Consider a system with n 
complex-valued amplitudes We split these amplitudes into real and imaginary parts 

and '0^, respectively, which evolve according to 


dt 

djA 

dt 


siW, 


(32a) 

(326) 


where 5^ and are the real and imaginary parts of the source function 5^. For this system 
the energy and enstrophy are given by 

E = (33) 

k 

and 

^ = (34) 

k 

where k ranges over the wavenumbers of the n modes. The properties of the source terms 
that guarantee conservation of energy and enstrophy are 


T,'^lSl + i’iSl=0, (35a) 

k 

+ = (355) 

k 

hence, we see that a system of n complex modes is completely equivalent to a system of 2n 
real modes. Therefore, the complex version of our conservative algorithms follows by applying 
the real algorithm separately to each component of the complex amplitudes. 
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D. Discussion 


It is worth saying a few words about computational efficiency. There are two sources of com¬ 
putational overhead associated with the conservative algorithms compared to the conventional 
methods. Here we concentrate on C-PC since C-Euler is not appropriate for practical use. 
In terms of operations, C-PC requires two additional multiplications and a square-root eval¬ 
uation over the standard predictor-corrector method. Importantly C-PC uses no additional 
storage. The cost of the extra operations will, in most cases, be negligible compared to 
the cost of one evaluation of Sf,. The square root is a cause for some concern as it may 
involve a function call. In any event, on modern hardware the square-root operation is only 
a small numbers of times (typically hve to seven) slower than multiplication. Furthermore, 
it is reasonable to expect that a conservative integrator will obtain a given global accuracy 
with a larger time step than the corresponding conventional integrator, thereby ameliorating 
the overhead problem. The second source of overhead is the occasional need to reduce the 
time step when the argument of the square root becomes negative. In practice we hnd that 
this happens approximately 10% of the time and, in light of the above discussion, we feel that 
this is not signfficant. 

Since the C-PC method reduces to the usual predictor-corrector algorithm in the small 
time-step limit, we expect that C-PC will inherit some of the stability properties of this 
method. While this does not establish the stability properties for large time steps, we hnd in 
practice that the C-PC method is numerically stable. 

In numerical studies of the Euler huid equations, an artihcial viscosity is often added to 
the dynamical equations to compensate partially for the spurious growth of the energy and 
enstrophy introduced by the numerical scheme. The viscosity is usually taken to vary as a 
power of wavenumber. However, only one of the two invariants can be exactly conserved 
by such a procedure and even this requires that the prescribed viscosity coefficient be time 
dependent. Moreover, this remedy can be shown to contaminate the modal evolution. In 
contrast, the conservative algorithms developed in this work faithfully reproduce the modal 
dynamics. 

In addition these methods can be applied to dissipative systems where the change in 
energy has a specihc physical origin. The same numerical errors that previously led to 
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nonconservation of energy will now contribute to the net energy change, thus having the 
effect of altering the underlying physics. For example, in a viscous fluid simulation the 
amount of energy leaving a mode is determined by the viscosity. It is a straightforward 
matter to use our methods to guarantee that the modal energy evolution is precisely that 
given by the physics. It is an open question and a subject of further investigation by the 
authors as to whether errors of this sort have the same structural effect on the solution as in 
the strictly conservative case. 

There is a simple interpretation of both the C-Euler and C-PC algorithms that sheds 
light both on their form and on the existence of the two branches, labeled by As numerous 
authors have observed, most traditional numerical methods conserve the linear invariants of a 
system. Consequently, one might be led to consider the possibility of transforming '0^ to new 
variables, in terms of which the invariants are linear. For the three-wave problem, this can 
be accomplished by making the transformation 4>fi = Upon applying the Euler method in 
the space and transforming back by taking the square root, one immediately obtains (16). 
This indicates that our restriction of the general constraint (13) to the condition (14) merely 
ensures that the modal energies evolve in a manner consistent with the Euler discretization 
of the energy equations. Below we will give derivations of integrators for other systems based 
on this idea. The C-PC algorithm can be viewed in the same light, except that the predictor 
is taken to have the simpler, nonconserving form. This also explains the '0^ = 0 fixed point 
in C-Euler: the modal energies have a second-order zero at '0^ = 0, thus it is no wonder that 
a first-order method fails at that point. 

IV. Lotka—Volterra 

As a further demonstration, consider the Lotka-Volterra predator-prey equations: 

^ = -jjx{l - y) , (36a) 

^ = y{l-x). {36b) 

These equations are surprisingly hard to integrate numerically since they are very susceptible 
to round-off error. With the exception of Kahan’s^^^l nontraditional method (which Sanz- 
Sernat^^l has shown to be symplectic) there are virtually no other methods that can integrate 
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this system without eventually failing due to round-off error. 

This is a noncanonical Hamiltonian systemt^^l with Hamiltonian 


and Poisson bracket 


H = X — logx -\- ny — julogy 

r, 1 (SS dg dg df 

{/.«} = xy 


( 37 ) 


(38) 


dx dy dx dy 

Jnst as with the three-wave problem, conventional integrators snch as Enler and predic¬ 
tor-corrector fail to conserve total energy H. It happens that the dynamics of this system 
are particnlarly sensitive to the valne of the energy, which explains the difficnlty that these 
methods enconnter when integrating (36). 

It is possible to derive a conservative algorithm for this systems nsing the methods of 
backwards error analysis ontlined above. The transcendental natnre of the fnnctions in the 
energy greatly complicates the procednre and prevents an analytical solntion of the rele¬ 
vant eqnations. In light of these problems, we take an alternative approach to deriving a 
conservative integrator. We proceed nsing the observation that standard methods snch as 
predictor-corrector exactly preserve linear invariants of a system of differential eqnations. 
To exploit this behavior, we introdnce new variables and ^2 defined by 


= a: — log a:, (39a) 

^2 = i^iy■ (396) 

This transformation was chosen so that H is a linear fnnction of and ^ 2 - Using the original 
eqnations of motion, we obtain 

^ = fi{x - l){y - 1), (40a) 

^ = -li{x -l){y-l). (406) 

Applying the nsnal second-order predictor-corrector to these eqnations yields 

'^i = ^i + Tiiix-l){y-l), (41a) 

^2 = ^i-r/i(a:-1)(|/-1), (416) 

^i(t + r) =^^ + ^fi[{x-l){y-l) + {x-l){y-l)] , (41c) 

^2(^ + ^) + • (41c?) 


20 




Strictly speaking, here x and y are to be computed from and .^2 by inverting (39a) and (396) 
respectively. Following the philosophy of the previous section, we instead compute x and y 
from the original equations of motion to obtain the following conservative integrator: 


X = X — T nx(l — y ), (42a) 

y = y + ry{l-x), (426) 

^i(t + r) = + - l){y - 1) + {x - l){y - 1)] , (42c) 

^2(^ + ^) =^2-^/^[(^-l)(?/-l) + (^-l)(^-l)] • (42c?) 


Here x(t -\- r) and y{t -\- r) are determined from + r) and + '^) by inverting (39). 
Since this inversion requires solving a transcendental equation, in practice it will have to be 
carried out iteratively. Although the expressions for x{t + T) and y(t + T) can not be written 
in closed form, (42) is still an explicit scheme. 

Notice that and ^2 ^ot one-to-one functions of x and y; has a minimum value 
of 1 at a: = 1, while ^2 has a minimum of // at y = 1. These minimum values play the same 
role as the point '0^ = 0 in the three-wave problem. Fortunately, the remedy is similar also: 
if either or ^2 are pushed below their respective minima, this indicates that the time step 
is too large. Temporarily reducing the time step alleviates this problem. 

To illustrate the effectiveness of our conservative algorithm, we integrate (36) taking y = 
1.5 with an initial condition of a:(0) = 1.0 and y(0) = 0.4. In Figure 7 we show a comparison 
between the standard predictor-corrector and C-PC. The C-PC orbit exactly conserves en¬ 
ergy and forms a closed curve. The predictor-corrector orbit spirals outward; a consequence 
of its energy gain. 

We provide this example to illustrate the generality of our method; however, since there 
is not an explicit expression for the inverse of the transformation (39) and a symplectic algo¬ 
rithm is known,this method seems to be of little practical value due to the computational 
overhead of iteratively determining x and y from and ^ 2 - 
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Figure 7: Integration of the Lotka-Volterra problem using a standard second-order predic¬ 
tor-corrector and the C-PC algorithm each with 8 x 10® time steps of size 0.02. A point is 
plotted every 200 time steps. The solid line represents the energy surface containing the initial 
condition. The points obtained from C-PC all lie on this curve. The dramatic effect of the 1.2% 
energy gain by the standard algorithm is clearly visible. 





V. Kepler Problem 


As a final example we consider the problem of a single particle moving in a gravitational 
potentialJ^^’^^^ Let r be the position vector of the particle of mass m and (j){r), where r = |r|, 
be the gravitational potential. The equations of motion for this system are 

dr 
dt 
dv 


= V . 


dt 


= -V(j). 


(43a) 

(436) 


This is a conservative system with the Hamiltonian 

H = -mv"^ -\- (j){r) 
2 


(44) 


As with all central force problems the total angular momentum, L = mr x v, is conserved, 
confining the motion to the plane perpendicular to L. We exploit this feature by aligning our 
coordinate system with the z direction parallel to L and introducing polar coordinates (r, 9) 
in the plane perpendicular to L. 

In these coordinates the equations of motion become 

dr 


dt 


dv^ 

dt 

dt 


f' 1 

m 
I 


m?r^ 


mr^ 


(45a) 

(456) 

(45c) 


where £ is the magnitude of the angular momentum and the Hamiltonian can be written as 

(46) 




Unlike all other central force problems, the Kepler problem has an additional constant of 
motion known as the Runge-Lenz vector. 


A = V X L + (hr . 


(47) 


This invariant is somewhat unique in that it is not the consequence of a symmetry of the 
equations of motion but is a result of the potential being proportional to 1/r. Conservation 
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of the Runge-Lenz vector can be associated with the fact that the orientation of the bound 
orbits of this system is fixed. It turns out the these orbits are elliptical and oriented with the 
major axis in the direction of A. We say “associated” here since in central force problems 
with any other force law, the orientation of the bound orbits processes. Furthermore, the 
Runge-Lenz vector is in some sense redundant: it is not needed to integrate the equation of 
motion, as there are already enough constants of motion to render the problem integrable. 

We adopt the initial conditions r(0) = Tq and 0(0) = n^(0) = 0, so that the vector A is 
in the a:-direction. Writing the potential as </)(r) = —K/r, where W is a constant, we see that 
the magnitude of A is given by 

A = - K. (48) 

mr^ 

A. A Conservative Integrator for the Kepler 
Problem 

The Kepler problem is an interesting example to consider in a study of conservative integra¬ 
tors, not only as a preliminary to studying multi-body problems (which are of astronomical 
significance), but also because of the Runge-Lenz vector. While this vector is functionally 
dependent on the Hamiltonian and on the angular momentum, exact conservation of these 
invariants neither guarantees conservation of the Runge-Lenz vector nor prevents the com¬ 
puted orbits from exhibiting a spurious precession. Hence numerical conservation of the 
Runge-Lenz vector is as much a structural issue as is conservation of energy. 

We now illustrate a conservative predictor-corrector (C-PC) for integrating (45) that 
exactly conserves H and A. The predictor is conventional: 

r = r -\- TV^ , (49a) 

Iff \ 

V^ = v +T 2- K] ^ (49ft) 

mr^ \mr ) 

0 = 0 + T-^. (49c) 

mr^ 


To obtain the corrector equations, we transform (r, n^) to the new variables 



(50a) 
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(506) 


^2 


1 2 1 
- 

2 2 mr 


2 ’ 


SO that H = + ^ 2 - Expressed in these new variables, the Hamiltonian is linear and will be 

conserved by conventional integrators. The corrector is given by 

^i(t + r) = ^1 + A, (51a) 

^ 2 (^ + ^) = ^2 - ( 51 ^) 


where 

A = ^ I 

2 I J’2 ^ 

In terms of the original variables, (51) may be rewritten as 



r{t + t) 


-K 

-Kjr + A ’ 


v^{t + T) = sgn(r;^) 






(52) 


(53a) 

(536) 


We still need an eqnation for 6. Thns far, we have enforced the invariance of H, bnt 
not A. Since only one integration variable remains to be determined, the conservation of A 
enforces the following constraint on 6: 


A ( cos 6 — 


mr 


sin0 = —K V , 


(54) 


as is seen upon taking the •u-projection of (47). To avoid the complexities associated with 
multiply-branched solutions, the most efficient method for solving (54) appears to be Newton- 
Raphson iteration, using 9{t) for the initial estimate. Convergence is rapid; typically, only 
3 or 4 iterations are required. 

In Figures 8 and 9 we present our integration results for the conventional C-PC and PC 
algorithms, respectively, adopting the initial parameters r = 1, 0 = 0, £ = 1, W = 3/2, 

and m = 1. To allow an even comparison, a slightly larger time step size was chosen for 
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X 


Figure 8: Solution of the Kepler problem computed using the conventional pre¬ 
dictor-corrector. A total of 1313 fixed time steps of 0.08 were used. 



Figure 9: Solution of the Kepler problem computed using the conservative pre¬ 
dictor-corrector. A total of 1000 fixed time steps of 0.105 were used. 
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the C-PC run such that the amount of computer time needed to reach the final time was the 
same in both cases. The new algorithm dramatically outperforms the traditional integrator. 
The artificial precession of the trajectory exhibited by the predictor-corrector result does not 
occur in the C-PC solution, due to the explicit conservation of the Runge-Lenz vector. 

B. Discussion 

We have demonstrated an explicit conservative integrator for the Kepler problem that cap- 
tnres all of its important strnctnral featnres. LaBndde and Greenspant^^’^^1 have constrncted 
integrators for this problem that conserve both energy and angnlar momentnm bnt it is 
nnclear whether their method exhibits orbital precession. 

The extension of these ideas to mnlti-body problems is a complicated task. Althongh 
each of the components of the angnlar momentnm are constants of motion, they are not in 
involntion. In the simple Kepler problem we are able to avoid any difficnlties associated with 
this behavior becanse the motion is confined to a plane perpendicnlar to the direction of the 
angnlar momentnm. In the mnlti-body problem this is no longer the case, which significantly 
complicates matters. 

VI. Conclusions 

We have demonstrated a method, based on the ideas of backwards error analysis, for deriving 
explicit, exactly conservative integration algorithms. This method consists of modifying the 
dynamical equations in such a way that when a particular conventional integration algorithm 
is applied to the modified equations, one obtains a solution consistent with the original 
equations that exactly conserves a system’s invariants. This method will generally yield 
an explicit algorithm when an explicit conventional algorithm is chosen as the basis for 
the conservative scheme. We have seen that this method can be interpreted in terms of a 
transformation to a new set of variables in which the invariants in question are linear. This 
promises to be a general method for deriving conservative integrators. 

In Section III, we saw that for a system with quadratic invariants, conservative integrators 
can be developed that are simple and computationally efficient. The case of quadratic invari¬ 
ants is of particular interest. The invariants in Lie-Poisson systems are typically quadratic 
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Casimirs. Furthermore, for Hamiltonian systems with Lie group symmetry, a Lie-Poisson 
system is the natural result of reduction; thus, our methods are applicable to integrating the 
dynamics on the Poisson manifold of such systems. For the integration of canonical Hamilto¬ 
nian systems where the configuration space is a Lie group Simo, Lewis and co-workerst^®^^®] 
have developed a series of methods that are symplectic and conserve momentum. One could 
imagine a hybrid of these algorithms: a conservative integrator of the type discussed above 
for integrating the dynamics on the Poisson manifold coupled to the algorithms of Simo et 
al. for reconstruction the full phase space flow. 

In addition to the desirable physical aspects of exact energy conservation there is some 
evidencet^®] that such conservation leads to nonlinear numerical stability of the algorithm. 
Furthermore, any conservative integration method developed for general systems conld cer¬ 
tainly be applied to Hamiltonian systems, providing an interesting comparison with symplec¬ 
tic methods. For example, this might shed some light on the choice between preserving phase 
space strnctnre and exact conservation of energy. In fact, one conld envision nsing the local 
change in phase space volnme as a diagnostic of the performance of a conservative integrator. 
These ideas will be the snbject of a fntnre paper. 
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